function [err]=EPSC_IO()
% Analyze Field Potential time course

err = 0;
EPSC_IO2(getmainselection, 0);
try
   sf = getmainselection;
   if(sf > 0) 
      pflag = getplotflag;
      QueMessage('EPSC_IO analysis', 1); % clear the que
      for i = 1:length(sf)
         EPSC_IO2(sf(i), pflag);
      end;
   end;
catch
   watchoff;
   QueMessage('Error in EPSC_IO Analysis routine', 1);
   err = 1
end;


function [err] = EPSC_IO2(sf, plot_flag)

global VOLTAGE CURRENT DFILE ALLCH
global CONTROL
%

err = 0;

[DFILE, err] = analysis_setup(DFILE, sf);

if(err ~= 0)
   return;
end;

% abstract general information
QueMessage('Preparing for analysis');
protocol=deblank(lower(CONTROL(sf).protocol));
rate = DFILE.rate.*DFILE.nr_channel/1000;

[records,pts]=size(VOLTAGE);
% always compute the time of each stimulus, in seconds from the start of the data
a = DFILE.ztime;
TM = (a - a(1))/1000; % time at which data was collected... gathered from ztime...
IO = seq_parse(CONTROL(sf).sequence);
if(isempty(IO))
   QueMessage('EPSC_IO: requires stim sequence', 1);
   return;
end;
IO = IO{:}';
% Get analysis windows 
if(ischar(CONTROL(sf).stim_time))
   pspt = CONTROL(sf).stim_time;
   p = find( pspt~= ' ');
   pspt = pspt(p);
   stim_time=seq_parse(pspt); % get the stimulus array - usually 2 times....
   stim_time = stim_time{:};
else
   stim_time = CONTROL(sf).stim_time;
end;

% get the psp definition array. For field potentials, we have the following
% sets of times: (layer 2 DCN pf FP)
% It is assumed that the same response is repeated with similar latencies for each stimulus...

if(ischar(CONTROL(sf).psp_time))
%   pspt = CONTROL(sf).psp_time;
%   p = find(pspt~= ' ');
%   pspt = pspt(p);
%   psp_time=seq_parse(pspt); % get psp definition array : for 
   psp_time = str2num(CONTROL(sf).psp_time);
else
   psp_time = CONTROL(sf).psp_time;
end;

% Get sign information
if(ischar(CONTROL(sf).sign_list))
   pspt = CONTROL(sf).sign_list;
   %p = find( pspt~= ' ');
   %pspt = pspt(p);
   % sign_list = seq_parse(pspt);; % get the stimulus array - usually 2 times....
   sign_list = str2num(pspt); % sign_list{:};
else
   sign_list = CONTROL(sf).sign_list;
end;

ts1=number_arg(stim_time(1,1)); % get the baseline duration (most of time before stimulus)
tpb=number_arg(psp_time(1:2:end));
tpe=number_arg(psp_time(2:2:end));

if(length(stim_time) > 1)
   ts2=number_arg(stim_time(1,2)); % get the baseline duration (most of time before stimulus)
end;


% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
RATES = (DFILE.rate .* DFILE.nr_channel) / 1000; % array of sampling rates, convert to msec
% calculate measurement times
fp_base = floor((ts1./RATES - 0.5)*0.75);
for i = 1:length(RATES)
   fpb(i,:) = floor(tpb./RATES(i));
   fpe(i,:) = floor(tpe./RATES(i));
end;
if(length(stim_time) > 1)
   dt = ts2-ts1;
   tpb2 = tpb + dt;
   tpe2 = tpe+dt;
   for i = 1:length(RATES)
      fpb2(i,:) = floor(tpb2./RATES(i));
      fpe2(i,:) = floor(tpe2./RATES(i));
   end;
end;

% measure holding current and "rmp"
for i = 1:records
   hold_cur(i) = mean(CURRENT(i,1:fp_base(i))); % save in array for later usage
   rmp_volt(i) = mean(VOLTAGE(i,1:fp_base(i)));
end;
CONTROL(sf).iHold = mean(hold_cur);
CONTROL(sf).Rmp = mean(rmp_volt);

[ismo] = FP_artsupp(CURRENT, DFILE, sf);

% but first smooth the CURRENT out a bit
for i = 1:records
   fsamp = 1000/RATES(i); % get sampling frequency
   fco = 3000;		% cutoff frequency in Hz
   wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
   if(wco < 1) % if wco is > 1 then this is not a filter!
      [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
      ismo(i,:) = DigitalFilt(b, a, ismo(i,:)')'; % filter all the traces...
   else
      ismo(i,:) = ismo(i,:);
   end
end

stdn = zeros(records, 1);
for i = 1:records
   base = mean(ismo(i,1:fp_base(i)));
   stdn(i) = std(ismo(i,[1:fp_base(i)]));
   ismo(i,:) = ismo(i,:) - base;
end;
stdn = sqrt(sum(stdn.^2)/records);
QueMessage('EPSC_IO - Measuring min EPSC current for stim 1');

MAX_S = zeros(records, size(fpb,2));
MAX_S2 = MAX_S;
twinpulse = 0;
if(length(stim_time) > 1 & all(fpe2 <= size(ismo,2)))
   twinpulse = 1;
end;
for i = 1:records
   for j = 1:size(fpb,2)
      if(sign_list(j) < 0)
         [npk, npi] = min(ismo(i,fpb(i,j):fpe(i,j)));
         % get a mean around 3 points...
         npi = npi+fpb(i,j)-1;
         npk = mean(ismo(i,npi-1:npi+1));
         if(j == 1)
            [ppk, ppi] = max(ismo(i,(fpb(i,j)+npi-1):fpe(i,j))); % neg to next positive (N1)
         else
            [ppk, ppi] = max(ismo(i,fpb(i,j):(fpb(i,j)+npi-1))); % first pos to main neg. (N2)
         end;
         ppk = 0;
         MAX_S(i,j) = abs(npk-ppk);
      else
         MAX_S(i,j) = max(ismo(i,fpb(i,j):fpe(i,j))); % simply positive
      end;
      if(twinpulse)
         if(sign_list(j) < 0)
            [npk, npi] = min(ismo(i,fpb2(i,j):fpe2(i,j)));
            if(j == 1)
               [ppk, ppi] = max(ismo(i,(fpb2(i,j)+npi-1):fpe2(i,j))); % neg to next positive (N1)
            else
               [ppk, ppi] = max(ismo(i,fpb2(i,j):(fpb2(i,j)+npi-1))); % first pos to main neg. (N2)
            end;
            MAX_S2(i,j) = abs(npk-ppk);
         else
            MAX_S2(i,j) = max(ismo(i,fpb2(i,j):fpe2(i,j)));
         end;
      end;
      
      
   end;
end;
%   MAX_P3S(i) = max(diff(vsmo(i,fp_p3b(i):(fp_p3b(i)+pk)'))) / RATES(i); % max slope on rising phase
%end;
for i = 1:records
   a = find(MAX_S(i,:) < stdn);
   allsignal{i} = a;
end;

maxiter = 500;
maxio = max(IO);
minio = min(IO);
if(exist('boltzanalysis'))
   for j = 1:size(fpb,2)
      mxs = max(MAX_S(:,j));
      sl = (mxs-min(MAX_S(:,j)))/maxio;
      fit_bv(j,:) = [0:maxio];
      [fp2, chisq2, niter2, fit_pki2, err2, dep2] = mrqfit('sum_boltzmann', ...
         [min(MAX_S(:,j)) mxs*0.5 sl maxio/5 mxs*0.5 sl maxio*0.5], IO, MAX_S(:,j), [], ...
         [0 1 1 1 1 1 1], [-1, -mxs*0.5, 0.0, -maxio*0.5, 0, 0, minio*2], ...
         [-1, mxs*2, 10, maxio*1.5, mxs*2, 10, maxio*1.5], maxiter);
      fit_ba_double(j,:) = sum_boltzmann(fit_bv(j,:), fp2);
      fit_ba_dbl(j,:) = fp2;
      if(fp2(2) < fp2(5))
         fpx = fp2(1:4);
      else
         fpx = fp2([1,5,6,7]);
      end;
      [fb, chisqb, niterb, fit_b, errb, depb] = mrqfit('sum_boltzmann', [min(MAX_S(:,j)) mxs sl*2 maxio], ...
         IO, MAX_S(:,j), ...
         [], [1 1 1 1], [-1, 0, 0.0, -maxio*0.5], [1, mxs*2, 10, maxio*1.5], maxiter);
      fit_ba(j,:) = sum_boltzmann(fit_bv(j,:), fb);
      f1{j} = fb;
      f1d{j} = fp2;
      
      
   end;
   
   if(twinpulse)
      for j = 1:size(fpb,2)
         mxs2 = max(MAX_S2(:,j));
         sl = (mxs2-min(MAX_S2(:,j)))/maxio;
         fit_bv2(j,:) = [0:maxio];
         
         [fp2_double, chisq2, niter2, fit_pki2, err2, dep2] = mrqfit('sum_boltzmann', ...
            fit_ba_dbl(j,:), IO, MAX_S2(:,j), [], ...
            [1 1 1 1 1 1 1], [-1, -mxs2*0.5, 0.0, -maxio*0.5, 0, 0, minio*2], ...
            [1, mxs2*2, 10, maxio, mxs2*2, 10, maxio*1.5], maxiter);
         f2d{j} = fp2_double;
         fit_ba_double2(j,:)=sum_boltzmann(fit_bv2(j,:), fp2_double);
         
         if(fp2_double(2) < fp2_double(5))
            fpx = fp2_double(1:4);
         else
            fpx = fp2_double([1,5,6,7]);
         end;
         [fb, chisqb, niterb, fit_b, errb, depb] = mrqfit('sum_boltzmann', [min(MAX_S2(:,j)) mxs2 sl*2 maxio], IO, MAX_S2(:,j), ...
            [], [1 1 1 1], [-1, 0, 0.0, -maxio*0.5], [1, mxs2*2, 10, maxio*1.5], maxiter);
         fit_ba2(j,:) = sum_boltzmann(fit_bv2(j,:), fb);
         f2{j} = fb;
      end;
   end;
   
else
   fit_bv=0;
   fit_ba=zeros(5,1);
   fit_bv2 = 0;
   fit_ba2 = zeros(5,1);
   f1 = 0;
   f2 = 0;
   f1d = 0; f2d = 0;
   fit_ba_double = zeros(maxio+1,1); fit_ba_double2 = zeros(maxio+1,1);
end;

EPSCIO.S1=MAX_S;
EPSCIO.S2=MAX_S2;
EPSCIO.IO=IO;
EPSCIO.S1_fitx=fit_bv;
EPSCIO.S1_fity=fit_ba;
EPSCIO.S1_fity2=fit_ba_double;
EPSCIO.S2_fitx=fit_bv2;
EPSCIO.S2_fity=fit_ba2;
EPSCIO.S2_fity2=fit_ba_double2;
EPSCIO.S1_par=f1;
EPSCIO.S1_par2=f1d;
EPSCIO.S2_par=f2;
EPSCIO.S2_par2=f2d;
EPSCIO.stdn = stdn;


CONTROL(sf).EPSCIO = EPSCIO; % copy the structure over

QueMessage('EPSC_IO - analysis complete');

na = size(MAX_S,2);
for i = 1:na
   IOX(:,i) = IO;
end;

% plot if figure is set
h = findobj('Tag', 'EPSC_IO'); % check for pre-existing window
if(isempty(h)) % if none, make one
   h = figure('Tag', 'EPSC_IO',...
      'Name', 'EPSC IO Analysis', ...
      'Units', 'normalized', ...
      'MenuBar', 'none', ...
      'WindowButtonMotionFcn', 'datac(''mouse_motion'', gcbf);', ...
      'WindowButtonDownFcn', 'datac(''mouse_down'', gcbf);', ...
      'WindowButtonUpFcn', 'datac(''mouse_up'', gcbf);', ...
      'NumberTitle', 'off');
end
figure(h); % otherwise, select it
clf; % always clear the window...
%Command Menu
uimenu('Label', 'Close &Window', 'Position', 1, 'callback', 'close(findobj(''Tag'', ''EPSC_IO''));' );
uimenu('Label', '&Print', 'Callback', 'print;');
uimenu('Label', 'Print&Setup', 'Callback', 'printdlg;');
fsize = 7;
msize = 3;
msg_pos = [0.85 0.90 0.15 0.07];
tmax = max(time);    

lab = {'1', '2', '3', '4'};
npm = size(MAX_S, 2); % number of measurments
spc = 0.03;
wid = 1/npm - spc;
if(npm == 1)
   wid = 0.5;
   spc = 0.15;
end;
for i = 1:npm
   subplot('Position',[spc+(i-1)*(wid+spc),0.550,wid,0.30]);
   hold on;
   P = plot(IOX, MAX_S(:,i), 's-', 'Markersize', msize);
   plot(fit_bv, fit_ba(i,:));
   plot(fit_bv, fit_ba_double(i,:), 'color', 'cyan');
   if(twinpulse)
      plot(IOX, MAX_S2(:,i), '+', 'Markersize', msize);
      plot(fit_bv, fit_ba2(i,:));
      plot(fit_bv, fit_ba_double2(i,:), 'color', 'magenta');
   end;
   
   set(gca, 'FontSize', fsize);
   ylabel('EPSC Amplitude (pA)');
   xlabel('Stim (uA)');
   grid;
   %legend(P, char(lab(i)));
   u=get(gca, 'YLim');
   ux = get(gca, 'XLim');
   datac_setcrosshair(gca, 'EPSCIO_info', 'uA', 'pA', msg_pos);
   
end;

 % plot the paired-pulse potentiation

if(twinpulse)
   PP = zeros(records, size(MAX_S, 2));
   for i = 1:records
      for j = 1:size(MAX_S, 2)
         if(MAX_S(i,j) < stdn)
            PP(i,j) = NaN;
         else
            PP(i,j) = MAX_S2(i,j)/MAX_S(i,j);
         end;
      end;
   end;
   subplot('Position', [0.55, 0.05, 0.43, 0.40]);
   plot(IOX, PP, 'o', 'Markersize', msize);
   hold on;
   %plot(fit_bv', fit_ba2'./fit_ba');
   plot([0:maxio], fit_ba_double2'./fit_ba_double');
   u=get(gca, 'YLim');
   set(gca, 'YLim', [0.5,u(2)]);
   grid;
   datac_setcrosshair(gca, 'EPSCIO_info', 'uA', 'Fac', msg_pos);
end;

% plot the raw data
subplot('Position',[spc,0.050,wid,0.40]);
%   plot(IO(1:records),MAX_P3S(1:records), '-r^', 'Markersize', msize, 'Markerfacecolor', 'r') % data in red
plot(time', ismo');
set(gca, 'FontSize', fsize);
v=get(gca, 'YLim');
vmax=max(abs(v));
%  set(gca, 'Ylim', [-vmax vmax]);
grid;
ylabel('(pA)');
datac_setcrosshair(gca, 'EPSCIO_info', 'ms', 'pA', msg_pos);


% plot textual information abstracted from analysis...
subplot('Position',[0.07,0.85,0.73,0.13])
r_uncomp = CONTROL(sf).access;
axis([0,1,0,1])
axis('off')
text(0,0.75,sprintf('%-12s R[%d:%d] %-8s',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol), 'Fontsize', 8);
text(0,0.5,sprintf('Solution:%-12s  gain:%4.1f  LPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'FontSize', 7);
text(0,0.25,sprintf('Noise: %8.3f', stdn));

orient landscape

% control printing and closing of window for automatic runs
% f = 1 creates plot and leaves it up
% f = 2 creates plot but closes it when done
% f = 3 creates plot and prints it and then closes it
% control printing and closing of window for automatic runs
% plot_flag = 0 creates plot and leaves it up
% plot_flag = 1 creates plot and prints to the printer
% plot_flag = 2 creates plot and prints it and then closes it
% checking the pflag box in datac main window is plot_flag = 1;
if (plot_flag > 0)
   h = findobj('Tag', 'EPSC_IO');
   hg = findobj('Tag', 'EPSC_info');
   set(hg, 'Visible', 'off');
   print(h);
   set(hg, 'Visible', 'on');
end
if plot_flag == 2  
   close
end

